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We investigate the non-equilibrium dynamics of the bosonic Hubbard model starting from inho¬ 
mogeneous superfluid or Mott insulator initial states using the truncated Wigner approximation 
(TWA). We find that the relaxation of the system develops in two steps for sufficiently large inter¬ 
action strengths: after a fast relaxation the system gets caught in metastable prethermalized states 
that precede the true equilibrium state. We find that the lifetime of these prethermalized states 
increases by several orders of magnitude as we increase the on-site interaction strength beyond a 
threshold value. We show that the emergence of long-lived metastable states in the quantum dy¬ 
namics is associated with an ergodic (active) to non-ergodic (inactive) dynamical phase transition 
in the ensemble of classical trajectories that contribute to the semiclassical limit. This dynamic 
phase transition, which is very similar to that found in different classical models of glasses, is closely 
related to the dynamic heterogeneity of the classical relaxation. 


I. INTRODUCTION 

The non-equilibrium dynamics of closed quantum sys¬ 
tems is nowadays a very active field of research. The mo¬ 
tivations are manifold. On one hand, the determination 
of the necessary conditions for a closed quantum systems 
to reach asymptotically a state of thermal equilibrium 
constitutes a key issue in the study of the foundations of 
statistical mechanics. On the other hand, with the ad¬ 
vent of cold atomic systems the manipulation of nearly 
perfectly isolated quantum systems with large coherence 
times has become experimentally feasible. There are sev¬ 
eral theoretical and experimental examples of systems 
that starting from non-equilibrium conditions reach, on 
accessible timescales, a state that is compatible with 
thermal equilibrium [1-4], at least at the level of simple 
correlation functions. However, situations in which ther¬ 
mal equilibration is elusive are clearly of particular inter¬ 
est. For integrable systems, thermal equilibration is not 
expected since the large number of constants of motion 
prevent the system from erasing memory of the initial 
conditions. However, there are examples of systems that 
are not integrable (but are close to an integrable point) 
for which thermal equilibration is not observed on the 
accessible timescales. Maybe the first experimental ex¬ 
ample of a (non-integrable) closed quantum many-body 
system failing to reach thermal equilibrium on accessi¬ 
ble timescales was provided by the famous experiment of 
Kinoshita et al. [5]. 

The phenomenon of prethermalization [6-13] is one of 
the possible explanations for the apparent lack of ther- 
malization in certain nearly integrable systems. Indeed, 
in certain situations it can be shown that the relaxation 
of systems close to an integrable point may proceed in 
two steps. The first step is a fast relaxation whose main 
mechanism is dephasing between some quasifree modes 
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of the system [10]. Close to an integrable point, where 
other relaxation mechanism such as inelastic collisions 
are extremely inefficient, such dephasing relaxation often 
gives rise to a highly non-thermal metastablc state. Af¬ 
terwards, the remaining relaxation channels become rele¬ 
vant and the prethermal metastable state starts to decay 
to the final state of the system [7, 8, 14]. One of the most 
characteristic phenomenons associated with prethermal¬ 
ization is the fact that certain observables reach its final, 
thermal value in the short timescales of the dephasing 
relaxation. The kinetic energy of the system is the most 
prominent example of such behavior. Prethermalization 
has been observed in experiments with ultracold bosonic 
gases [11, 15]. 

In particular, prethermalization may be behind the ap¬ 
parent breakdown of thermalization in the strongly inter¬ 
acting Bose-Hubbard model, which has received special 
attention. Starting with the pioneering work of Kollath 
et al. [16], where using density matrix renormalization 
group (DMRG) techniques it was shown that for suffi¬ 
ciently strong interactions few-body correlation functions 
reached a quasistationary state that was not compatible 
with thermal equilibrium, several studies have focused 
on this system [17, 18]. In Ref. [18] it was shown that 
for strong interactions the system gets trapped in long- 
lived inhomogeneous metastable states. Such kinetic ar¬ 
rest, was shown [18] to be associated with a dynamical 
localization in the many-body Hilbert space. The Bose- 
Hubbard model is also one of the most extensively real¬ 
ized models in cold atoms experiments [19]. In particular, 
experiments with inhomogeneous Bose-Einstein conden¬ 
sates also demonstrate pronounced timescale separation 
and slow thermalization for strong interactions [20-22]. 

On the other hand, since the seminal work of Sred- 
nicki about the eigenstate thermalization hypothesis [23] , 
semiclassical considerations have been very important in 
our understanding of quantum thermalization. Recently 
Cosme et. al [24] made a contribution in this line through 
the truncated Wigner approximation (TWA) applied to a 
two sites multi-band bosonic Hubbard model. Since the 
TWA allows to consider the first order quantum correc- 
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tions to the classical action [25] the connection between 
quantum physics and its classical limit becomes explicit. 
In particular, it has been checked in [24] that the ther- 
malization of the quantum system comes along the er- 
godicity of classical trajectories. Also along these lines, 
of particular relevance to the study presented here are 
the works [26] and [27] where the breakdown of thermal- 
ization for strong interactions was linked with the chaotic 
properties of the mean-field limit. More precisely, it was 
shown that the Lyapunov exponent of the classical tra¬ 
jectories associated with the solutions of the mean-field 
equations is suppressed for strong interactions. 

The appearance of long-lived prethermalized states 
in the relaxation of isolated quantum systems is phe¬ 
nomenologically very similar to the physics of classical 
glassy systems such as spin glasses [28], atomistic glass 
formers [29, 30] and kinetically constrained models [31]. 
Such systems exhibits a typical two-step relaxation: after 
a fast inertial regime (analogous to the dephasing regime) 
they get caught in metastable states whose lifetime can 
be astronomical in certain parameter regions, for exam¬ 
ple, for sufficiently low temperature or high density in 
the case of atomic models. One of the most distinctive 
characteristics of glassy relaxation is the appearance of 
dynamical heterogeneity [31], i.e., the fact that the re¬ 
laxation is spatially non-homogeneous: fast and slow re¬ 
gions are clustered together. One of the most appealing 
theories to explain such feature is that of a dynamical 
phase transition between active and inactive phases tak¬ 
ing place in the ensemble of possible histories that the 
system can follow in the relaxation [30, 32-34]. 

In this work we revisit the problem of the appar¬ 
ent breakdown of thermalization in the bosonic Hub¬ 
bard model from the semiclassical perspective provided 
by the TWA. For large occupation numbers, this ap¬ 
proach allows to access the dynamics of the quantum 
system performing weighted averages over ensembles of 
classical (mean field) trajectories. We explicitly show 
that prethermalized states can have a huge lifetime 
that grows exponentially while increasing the coupling 
strength, which prevents the observation of thermal equi¬ 
libration on any reasonable timescale. Moreover, we find 
that the emergence of long-lived prethermalized states is 
closely related with the appearance of non-ergodic, glassy 
features in the mean-field trajectories associated with the 
quantum dynamics. More specifically, we find that the 
ensemble of mean-field trajectories undergo a dynamical 
phase transition between an active and an inactive phase, 
completely analogous to that observed in classical models 
of glasses. We find that in the present case this transi¬ 
tion is also intimately correlated with a remarkable phe¬ 
nomenon present in the mean field trajectories, dynam¬ 
ical heterogeneity. We shall go deeper into these issues 
in the remaining part of the paper, which is organized 
as follows: In Section II we describe the model, the ini¬ 
tial conditions, the semiclassical approximation used to 
study the quantum dynamics and show the results for the 
quantum dynamics, while in Section III we thoroughly in¬ 


vestigate the properties of the ensemble of mean-field tra¬ 
jectories that contribute to the quantum dynamics and 
show that they show clear signatures of glassiness as we 
increase the strength of the interaction. In Section IV we 
briefly sumarize our results. 


II. QUANTUM DYNAMICS 

Let us consider the ID bosonic Hubbard model with 
Hamiltonian H = Hq + Hi nt : 

L 

Hq = ( a \ a j +1 + a \+l a j) ) 

3 = 1 
L 

Hint — U ^ > (-0 

j=1 

where L is the number of sites in the chain, the a’s are 
canonical bosonic operators and rij = a'-a :) . We shall 

consider periodic boundary conditions, a# +m = a#. We 
shall analyze the dynamics generated by (1) starting from 
an inhomogeneous initial state described by a density 
matrix po , which we left unspecified until the next sec¬ 
tion. In particular, we shall focus on the relaxation of 
the density profile {nj{t)) = Tr [p 0 e' lHt nje~ lHt ], where 
we have set K = 1. The thermal density profile is given 
by {nj)th = N a /L = N for all j, where N a is the total 
number of bosons. 

When the average number of atoms per site N is large 
the dynamics of the Bose-Hubbard model can be studied 
using the TWA [35] . To calculate the time dependence of 
expectation values within the TWA we have to consider 
the solution of the classical equations of motion associ¬ 
ated with the quantum Hamiltonian (1). These are the 
standard lattice Gross-Pitaevskii equations [35], 


i -= -J&j+iit) +Vb'-iW) + 2 Uipj(t) \^j{t)\ 2 .(2) 

where the classical fields are normalized to the total num¬ 
ber of particles IVb'WI 2 = N a . Then the expec¬ 

tation value of any given operator H at time t can be 
calculated averaging the corresponding classical observ¬ 
able fl c i over an ensemble of initial conditions weighted 
according to the Wigner transform of the initial density 
matrix po 


Q(t) = / #S#oP(V’o,V’o) , (3) 


where p is defined as 

P(i> o.V’o) = J drjodrio (i>o~ Po V>o + y) 

x e -\^o\ 2 -z\vo\ 2 ^(Vo^o-Vo^o) _ 


( 4 ) 





3 


To lighten the notation we have omitted the site index j. 
The measure is 

L 

#o#o = II difjj (0(0). (5) 

3 = i 

The correspondence between classical and quantum 
observables can be formulated most easily using the co¬ 
herent state Bopp representation, which makes the as¬ 
signments 


, /* 19 

a ^ * ~2ty' 

(6) 

1 d 


“ -*■ ' P+ 2W 

(7) 

"j ->■ IVbl 2 - l ■ 

(8) 


The TWA is the leading order approximation in an ex¬ 
pansion in a small parameter that measures deviations 
from classicality [25]. In the particular case of interact¬ 
ing bosons this parameter is 1/-/V, the inverse of average 
number of bosons per site. The TWA is a controlled ap¬ 
proximation in the sense that it is possible to calculate 
higher order corrections that take the form of stochastic 
perturbations to the classical trajectories [25]. However, 
as any approximation, be it controlled or not, there are 
certain parameter regimes in which it works better than 
others. To define when we expect this approximation 
to be accurate we may introduce the nonlinearity pa¬ 
rameter A = that is the ratio between the typical 
potential energy per site and the typical kinetic energy 
per site [27]. When the interactions are strong enough, 
A ~ N 2 , the system (in equilibrium) undergoes a quan¬ 
tum phase transition to a Mott insulating state. In the 
vicinity of the transition quantum fluctuations become 
large and we cannot expect the TWA to work well in this 
case, i.e., second order corrections become important in 
that regime. However, there is a wide regime A < N 2 
(U / J < N) where the ground state of the system is a 
superfluid (weakly or strongly interacting depending on 
whether A < 1 or A > 1 respectively) where we can expect 
the TWA to be a good approximation to the dynamics. 
We shall work in that regime. 

In parallel, the choice of the initial condition is also 
relevant for the accuracy of the results. We shall work 
with two types of initial states, coherent (superfluid) ini¬ 
tial states and Fock (Mott insulator) initial states. The 
coherent state represents the initial condition that is clos¬ 
est to a perfectly defined classical initial condition and 
therefore we can expect the TWA to accurately capture 
the dynamics. In particular, the Wigner function is 

Pc(^r 0 ) = U 2 -^- Nj ^ ( 0 ) 

j. j. 7 t 

j= 1 


where Nj = (n,j( 0)). The Wigner function is a true 
probability distribution and can be sampled efficiently 
as Vb = \fN~o + \ + * 772 ) [24, 36] where ?y 1 and rj 2 are 

two real Gaussian random variables. 




FIG. 1: Dynamical distance d(t ) as a function of time for the 
Fock (up) and coherent (down) initial states. Different curves 
correspond to different interaction strengths U/J. 

The Fock initial state is the least classical initial state 
in the sense that the initial phases are completely un¬ 
specified. The Wigner function in this case is 

L 

Pf(M S) = II 2e- 2| ^ |2 T iV . (4 |^| 2 ) , (10) 

j =1 

where Ln(x) is the N th order Laguerre polynomial. Un¬ 
fortunately, this p generically is not definite positive and 
for large N has a highly oscillatory behavior. This fact 
makes convergence much slower when doing a Monte 
Carlo integration over initial conditions. To bypass this 
issue we shall approximate this Wigner function by prod¬ 
ucts of true distributions reproducing the first moments 
of the Laguerre polynomials. The simplest way to do 
this is just to reproduce the first moment by replac¬ 
ing each factor ffv, (ip*. ijjj ) by a Dirac delta function 
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S (j'lpj | 2 — \ — Nj) and averaging over random phases 
when doing the Monte Carlo integration [24, 36, 37]. The 
next step would be considering a distribution with two 

moments, i.e., the Gaussian —. We 
will refer to this two approximations as the delta and 
Gaussian Wigner functions p$ and p g respectively. We 
have checked that in the particular cases that we inves¬ 
tigated both approaches give very similar results. We 
shall therefore use the delta function approximation due 
to computational convenience. For the Fock initial state 
it was found that, in the worst case, the TWA results 
quantitatively deviate from available exact solutions for 
times larger than t c = J/U [35]. However, the predic¬ 
tions of the TWA remain qualitatively valid for larger 
times. 

We shall now proceed to investigate the quantum dy¬ 
namics of the density profile. We shall chose the coherent 
or Fock initial states in such a way that the initial density 
profile is (20, 0,20,0 ,...) : an alternation of empty and 
highly occupied sites. With this election N = 10. The 
system is on a ring of size L = 30. This initial state may 
be experimentally relevant for cold atomic gases loaded 
on optical lattices [1]. To calculate the quantum dynam¬ 
ics we sample 10 3 different initial conditions both for the 
coherent and for the Fock initial state, but in order to test 
convergence we went up to 10 4 realizations. To quantify 
the overall relaxation of the density profile to its thermal 
configuration we introduce the dynamical distance: 

d?(t) = - N f ■ (11) 

3 =1 

Clearly, thermal equilibrium implies d 2 (t) = 0. In Fig. 1 
we show the decay of d 2 (t) for both the coherent and 
the Fock initial states for different values of the coupling 
strength. In both cases for low coupling strength the 
system thermalizes quickly in a timescale of the order of 
one hopping [1, 2], As we increase the coupling strength a 
metastable state emerges in between the initial relaxation 
and the final decay to equilibrium. The lifetime of such 
metastable state increases with the coupling strength un¬ 
til it becomes larger than the maximum time available in 
our simulation. For the Fock initial state the increase 
of the lifetime of the metastable state is more abrupt. 
In both cases, the density profile of the metastable state 
is closer to that of the initial state as we increase the 
coupling strength. 

In Fig. 2 we show the relaxation of the kinetic en¬ 
ergy ek in (t) = Ti[p 0 e lHt H 0 e~ lHt ] for the coherent state 
initial condition for several values of the interaction 
strength. For the Fock initial state we obtain qualita¬ 
tively the same results. We observe that for weak in¬ 
teractions U/J < 0.1 the kinetic energy prethermalizes 
in a timescale tp/ ak ~ 10J -1 . For intermediate interac¬ 
tion strengths 0.2 < U/J < 0.7 the kinetic energy also 
tends to a well defined constant but with a larger relax¬ 
ation time of the order of ~ 30J -1 ( int stands for 
intermediate). Whereas for strong interaction strengths 



FIG. 2: Relaxation of the kinetic energy for the coherent 
state initial condition. Different curves correspond to different 
interaction strengths U/J. Curves are vertically shifted for 
clarity. 


U/J> 0.7 the prethermalization timescale is t *‘ rons ~ 
3 J^ 1 . The variety of prethermalization timescales can 
be understood from the fact that the quasifree modes 
behind the dephasing relaxation are qualitatively differ¬ 
ent for strong and weak interactions. For weak interac¬ 
tions the quasifree modes are related to the momentum 
eigenmodes that diagonalize Hq and are completely de¬ 
localized in space. For strong interactions, the quasifree 
inodes are related with the eigenmodes that diagonalize 
Hi n t and are completely localized excitations. For inter¬ 
mediate interactions there is a truly non-trivial regime 
for which the quasifree modes are neither completely lo¬ 
calized nor delocalized. This suggests that the effective 
models that describe the short time dynamics in the three 
cases are completely different. However, a remarkable 
fact is that for strong interactions the prethermalization 
time is the same regardless the specific value of the cou¬ 
pling strength, while the relaxation timescale of the den¬ 
sity profile shows large variations (see Fig. 3) due to the 
presence of metastable states. 

The crossover from fast thermalization to the 
metastable state dominated regime for the two different 
initial conditions occur in different coupling ranges. To 
obtain a decay timescale of d 2 (t) we fit an exponential 
to the tail of the decay. In Fig. 3 we show the results 
for the decay timescales for both initial conditions. We 
find that the decay timescales increase by several order 
of magnitude as we increase the coupling strength in a 
small range. We find that the center of the crossover 
region is around U/J ~ 0.7 for the Fock initial state 
and U/J ~ 1 for the coherent state. The emergence of 
prethermalized states with such large lifetime may com¬ 
pletely hinder the observation of the thermal equilibrium 
state of the system after the quench and certainly is at 
the root of the apparent lack of thermalization observed 
in earlier works [16, 38]. We find the situation rather 
similar to that of glasses, systems that exhibit the typi- 
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FIG. 3: Logarithm of the relaxation time as a function of 
coupling U/J for the Fock (up) and coherent (down) initial 
state. 


cal two-step relaxation as a consequence of getting caught 
in extremely long-lived metastable states, whose lifetime 
can be astronomical for sufficiently low temperatures or 
high densities. In the next section we shall outline a more 
qualitative relation between the emergence of long-lived 
metastable states and glassiness by analyzing the prop¬ 
erties of the mean field trajectories associated with the 
quantum dynamics. 


III. GLASSY PROPERTIES IN THE 
MEAN-FIELD TRAJECTORIES 

In Ref. [24] it was proposed that the ability of a quan¬ 
tum system to thermalize is related with the ergodicity 
of the classical trajectories of the associated mean-field 
system, an idea that is also implicit in previous studies, 
such as [26] and [27]. In this Section we shall show that 
the emergence of long-lived prethermalized states in the 
quantum dynamics of the Hubbard model is correlated 
with the lack of ergodicity of the classical trajectories 
that are used to build up the quantum dynamics. 

To qualitatively introduce the discussion, in Fig. 8 we 
show a stroboscopic sampling of the position in phase 
space of the field ipi for two different couplings. We can 
see that for small coupling {U = 0.5) the system explores 
the whole phase space while for large coupling (U = 1.5) 
it remains captured in a ring of finite width in the Re^i) 


FIG. 4: Position in the phase space Re^i) vs. Im(-i/>i) for 
the first site of the chain observed every a J At = 0.1 for 
J tobs = 40 for U/J = 0.5 (up) and U/J = 1.5 (down). The 
big red dot corresponds to the initial configuration. Each site 
was chosen to have a random phase as initial condition. 


vs. Im( , !/>i) space. A similar behavior is observed when 
analyzing the behavior of the phase space trajectory of 
the other sites. Going back to Fig. 1 we can check that in¬ 
deed for U = 0.5 the quantum density profile is thernral- 
ized for tJ = 40 while it is still trapped in the metastable 
state for U = 1.5. 

Additionally, we point out that the phase of the fields 
0j(t) = arctan(Im['i/>j(t)]/Re[^(t)]) turns out to be al¬ 
ways ergodic, in the sense that along the dynamic evo¬ 
lution it visits all values between 0 and 27r, irrespective 
of the value of the coupling strength, for all trajecto¬ 
ries. This fact is also illustrated by Fig. 4 where, even 
for the strong coupling regime, the phase of the field 
explores all values. The lack of ergodicity is thus en¬ 
coded in the dynamics of the particle density at each site 
nj(t) = | 2 . 

We will now proceed to make a more detailed and 
quantitative analysis of the properties of the classical tra¬ 
jectories. In order to characterize them we introduce the 
mobility K of each trajectory 

tabs n 

k m] = ^ E E (i^(*+ Ai )i 2 - (i^coi 2 ) 2 , (i2) 

t=0 j =1 

where t 0 f, s is the maximum observation time, i.e., the 
extension of the trajectories and At is a UV cutoff that 
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kills possible small oscillations leaving just the actual dis¬ 
placement in the phase space. K is an extensive- 

in-time quantity that will typically be large for ergodic 
trajectories, while it will be generically small for non- 
ergodic trajectories as can be intuitively seen from Fig. 
8 . Thus, K [ijj(t)] can be used as an order parameter to 
distinguish ergodic (mobile) from non-ergodic (immobile) 
trajectories. 

Since we are interested in assessing the properties of 
the trajectories that are most relevant to the quantum 
motion we will introduce a measure on the space of clas¬ 
sical trajectories given by the Wigner function of the 
quantum initial state that we wish to consider. In the 
case of the coherent state, the measure pc(il> o> V’o) ' s we ^ 
defined while in the case of the Mott insulating state we 
can choose any well defined distribution that approxi¬ 
mates pf(V’o,'0o)> suc h as p$ or p g . For simplicity, we 
shall restrict in the following analysis to the coherent 
state case. 



0 100 200 
K/btobs-l 

FIG. 5: Histogram of the mobility K[ip(t)] for ensembles of 
2 x 10 3 trajectories sampled according to the Wigner func¬ 
tion of the coherent state initial condition pc(ipo, ipo) f° r f° ur 
different values of U/J. 


In Fig. 5 we show the mobility histogram from an 
ensemble of 2 x 10 3 trajectories sampled according to 
PC'W’o,V , o) f° r f° ur different values of U/J. The center 
of the distributions monotonically shifts to lower values 
of the mobility as we increase U/J , but the dispersion 
around the center decreases for large coupling constants. 
In particular, for sufficiently large coupling all the tra¬ 
jectories have a very small mobility with a small disper¬ 
sion, showing that all trajectories are freezed, at least 
up to timescales of the order of t 0 b s ■ However, there 
is an intermediate coupling regime 0.3 < U/J < 0.6, 
that corresponds to the intermediate regime described 
when analyzing the relaxation of the kinetic energy in 
the previous section, for which the dispersion almost 
doest not changes. This confirms that inactive, non- 
ergodic classical trajectories are related with the exis¬ 
tence of non-thermal metastable states in the quantum 
dynamics. Moreover, since the shift of the center of the 


histograms is continuous we may expect that for inter¬ 
mediate couplings there should be a coexistence of active 
and inactive trajectories, in which case the mobility of 
the trajectories would depend on the details of the ini¬ 
tial conditions. This is indeed the case. This observa¬ 
tion can be cast in the language of a dynamical phase 
transition [30, 32-34, 39] taking place in the ensemble of 
classical trajectories. To show this we first construct a 
canonical probability distribution in such ensemble, cou¬ 
pling the extensive-in-time order parameter K [il’(t)] with 
an intensive field s. This distribution will be proportional 
to 


P s [m]<xPo[m]e-' Kmt)] , (13) 

where the Pq [ijj(t)] is the s = 0 probability distribution 
that we take as the Wigner function of the coherent state 
Pci^o^o) as discussed earlier. 

We will compute expectation values in the ensemble 
of trajectories by summing over all trajectories weighted 
with P s [ip(t)] in the following fashion 

= (ci [v>(t)]) s = y [V’W] n bP(t )]» (14) 

s m 

where Cl is a trajectory functional and Z s is the 
partition function 

z s = j2 p sim}- (is) 

b(t) 

This way of averaging will give a different weight to mo¬ 
bile and immobile trajectories depending on the value of 
s: for larger s more relative weight is assigned to inactive 
trajectories and viceversa. 

In Fig. 6 we show the average order parameter K s as 
a function of s for U/J = 0.53, a value of U/J for which 
we have observed coexistence of active and inactive tra¬ 
jectories. The behavior of K s mimics the behavior of 
the order parameter in a finite volume equilibrium phase 
transition: it shows a marked step between two well de¬ 
fined values corresponding to the two different phases 
in the ensemble (active and inactive phases) and, more¬ 
over, while increasing t 0 b s (the analogous of the volume 
in equilibrium) the step in K s becomes more and more 
sharp. This can be appreciated more clearly looking at 
the susceptibility: 

r)K 

xs = —gf = mm-K s ) 2 )s, (i6) 

which exhibits a peak that grows while increasing t 0 bs- A 
similar scaling behavior can be observed with increasing 
L. The critical value is located around zero, s* = 0. This 
dynamical phase transition is very similar to that found 
first in kinetically constrained models of glasses [32] and 
then in atomistic models of glass formers [34] and quan¬ 
tum systems [39]. The main difference that we observe 
with respect to the works [32, 34, 39] is that the present 
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FIG. 6: Dynamical phase transition in trajectory space. Up¬ 
per panel: Averaged order parameter K s as a function of the 
intensive field s for an ensemble of 10 3 trajectories sampled ac¬ 
cording to the coherent state Wigner function for U/J = 0.53. 
Lower panel: Dynamical susceptibility y s . 


dynamic phase transition does not seems to be of first or¬ 
der. In a first order (dynamic or static) phase transition 
the distribution of the order parameter at coexistence is 
bimodal due to surface tension effects between the coex¬ 
isting domains of different phases. In the case analyzed 
in this work the order parameter distribution at coex¬ 
istence (s ~ 0 for U/J = 0.5) is unimodal, as can be 
inferred looking at Fig. 5, which is compatible with a 
continuous phase transition. 

To better characterize the physical properties of the 
phases we will discuss an ergodicity measure. We define 
the overlap correlation function, 

L 

Cs(t ) = £<(«,-(*) - A0K-(0) - N)) s . (17) 

j=1 

C s ( t ) quantifies the overlap between the configuration at 
time t and the initial density configuration. The extent 
to which it is non-zero in the limit of large t is a mea¬ 
sure of nonergodicity. In Fig. 7 we show the correlation 



FIG. 7: Correlation function C s (t) as a function of time for 
U/J = 0.53 in the active (s < 0) and inactive (s > 0) phases. 
We also show, for comparison, the correlation function at co¬ 
existence (s = 0). 


function C a ( t ) for the active (ergodic) and inactive (non- 
ergodic) phases. In the active phase (s < 0) the correla¬ 
tion function rapidly relaxes to a small value, while for 
the inactive phase (s > 0) trajectories remain trapped in 
a single state throughout the observation time. 
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FIG. 8: False color plot of the local mobility Qi(t) for two 
representative trajectories corresponding to U/J = 0.6 (upper 
panel) and U/J = 0.1 (lower panel). 


In classical models of glasses such dynamical phase 



























transition represents one of the main theories to account 
for an striking and distinctive feature of glassy materi¬ 
als: dynamic heterogeneity [30]. In contrast with normal 
fluids, the relaxation of glasses is heterogeneous in space: 
fast and slow regions are clustered together. We found 
that individual mean field trajectories (that correspond 
to single solutions of the Gross-Pitaevskii equation) ex¬ 
hibit the phenomenon of dynamical heterogeneity. We 
shall define a local measure of mobility Qi(t) as 

Qi(t ) = I rii(t) - m( 0)1, (18) 

that measures how much the in site density differs from 
the initial value. In Fig. 8 we show density plots of Qi(t) 
for representative trajectories corresponding to U/J = 
0.1 and 0.6. For U/J = 0.1 (the weak coupling regime) 
the dynamics of the system, as captured by the indicator 
Qi(t), is fairly homogeneous. For U/J = 0.6 (entering the 
strong coupling regime), sites with high local mobility 
and low local mobility are clustered together in space- 
time. This is the signature of dynamic heterogeneity. It 
is remarkable that the simple Gross-Pitaevskii equation 
is able to generate such complicated dynamics. 

The relation between the dynamical phase transition 
and the dynamical heterogeneity can be explained using 
a simple reasoning: the dynamical phase transition pic¬ 
ture implies that for intermediate values of U/J different 
initial conditions may trigger slow or fast dynamics, de¬ 
pending only on the details of the initial condition; if 
we consider a sufficiently large system with short range 
interactions, distant regions in space will not be corre¬ 
lated and, if one region has a “fast” initial condition and 
other a “slow” initial condition, dynamical heterogene¬ 
ity may arise. Moreover, by the same argument, during 
the dynamical process slow regions may turn into fast 
regions and viceversa. On contrary to what is observed 
in [40, 41], the quantum dynamics does not show dy¬ 
namical heterogeneity. The average needed to obtain the 


quantum dynamics from the individual classical trajecto¬ 
ries washes out any initial condition dependent dynam¬ 
ical heterogeneous pattern, which is associated with the 
fact that the density profile has a well defined relaxation 
timescale as seen before. 


IV. CONCLUSIONS 

We have studied the dynamics of the Bose-Hubbard 
model following a quantum quench from an inhomoge¬ 
neous initial state. Using the TWA we have analyzed 
the dynamics of the density profile of the system, that 
exhibits long-lived prethermalization plateaus for mod¬ 
erate and large couplings. We were able to relate this 
fact with the non-ergodic propoerties of the associated 
mean-field system. In particular, we have shown that 
the ensamble of mean field trajectories that are relevant 
to the quantum motion undergoes a dynamic phase tran¬ 
sition from an active to an inactive phase, much alike to 
what is found in classical models of glasses. We have also 
shown that this transition is related to the phenomenon 
of dynamic heterogeneity in the classical trajectories, one 
of the hallmarks of glassiness. In sum, we believe that 
our work opens a new perspective on the dynamics of 
closed quantum systems by relating the physics of clas¬ 
sical glasses to the well known phenomenon of prether¬ 
malization. 
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